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ABSTRACT 

A quarter of DA white dwarfs are metal polluted, yet elements heavier than he- 
lium sink down through the stellar atmosphere on timescales of days. Hence, these 
white dwarfs must be currently accreting material containing heavy elements. Here, 
we consider whether the scattering of comets or asteroids from an outer planetary sys- 
tem, following stellar mass loss on the asymptotic giant branch, can reproduce these 
observations. We use N-body simulations to investigate the effects of stellar mass loss 
on a simple system consisting of a planetesimal belt whose inner edge is truncated by 
a planet. Our simulations find that, starting with a planetesimal belt population fitted 
to the observed main sequence evolution, sufficient mass is scattered into the inner 
planetary system to explain the inferred heavy element accretion rates. This assumes 
that some fraction of the mass scattered into the inner planetary system ends up on 
star-grazing orbits, is tidally disrupted and accreted onto the white dwarf. The sim- 
ulations also reproduce the observed decrease in accretion rate with cooling age and 
predict accretion rates in old (>lGyr) white dwarfs, in line with observations. The 
efficiency we assumed for material scattered into the inner planetary system to end up 
on star-grazing orbits is based on a Solar-like planetary system, since the simulations 
show that a single planet is not sufficient. Although the correct level of accretion is 
reproduced, the simulations predict a higher fraction of accreting white dwarfs than 
observed. This could indicate that evolved planetary systems are less efficient at scat- 
tering bodies onto star-grazing orbits or that dynamical instabilities post-stellar mass 
loss cause rapid planetesimal belt depletion for a significant fraction of systems. 

Key words: planetary systems, Kuiper belt, white dwarfs, planets and satellites: 
dynamical evolution and stability, (stars:) circumstellar matter 



1 INTRODUCTION 

Keck observations of cool single DA white dwarfs find that 
~25 % contain elements he avier than helium in their spec- 
tra jZuckerman et al]|2003l '). These elements sink rapidly in 
the white dwarf's atmosphere and their presence means that 
these white dwarfs must be currently accreting material con- 
taining heavy elements. Initially it was thought that these 
observations were a signature of accretion from the interstel- 
lar medium; however, this was ruled out by a lack of corre- 
lation between their accreted calcium abundances and spa- 
tial kinematical distributions relative to interstellar material 
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comets from the remnants of main sequence planetary sys- 
tems are scattered onto orbits that approach close to the 
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star, due to altered dynamics following stellar mass loss on 
the asymptotic giant branch. Bodies that come within the 
tidal radius of the star are disrupted, potentially forming 
a dusty disc, before accreting onto the star. Spitzer obser- 
vations of some of the most highly polluted systems that 
find excess emission in the near-infra-red, consistent with a 
close-in dusty disc. Such a disc is observed aro und 1-3% of 
white dwarfs with cooling ages less than 0.5Gyr (jFarihi et al.l 
120091 ). 

Although the disruption of an asteroid or comet is 
widely quoted as the explanation for such systems, the fea- 
sibility of this process has not been thoroughly investigated. 
Evidence that the accreted material is asteroidal in nature 
is high; the co mposition of the a ccreted material in systems 
such as GD40 (|Klein et al.ll2010T ) highly resembles asteroids 
in our solar system. Firstly, in order for this to be the case, 
planetesimals must survive the star's evol ution . Consi dering 
only stellar wind drag and sublimation, lJural (|2008l ^ show 
that asteroids of l-10km i n size survive the giant branch 
evolution outside of 3-4AU. iBonsor fc WvattJ (|2010l ) model 
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in more detail the evolution of the observed population of 
debris discs around main sequence A stars, showing that 
white dwarfs should have planetesimal belts, but that these 
are hard to detect. 

A mechanism is still required to transport material from 
an outer planetesimal belt to the star. Although other sug- 
gestions, such as stellar encounters (Farihi et al. 2010 in 
prep), have been made, scattering by planets is the most 
likely mechanism. It is not clear whether planets will sur- 
vive the star's evolution. Many close-in planets will be swal- 
lowed by the expandi ng stellar envelope whil st the star is on 
the giant branch (e.g. lVillaver fc Livicl (|2007f )). Multi-planet 
systems m ay also become dynamicall y unstable post-stellar 
mass loss iDebes fc Sigurdssonl 120021 ). and planets may be 
ejected or collide with the star. 

The dynamics of multi-planet systems post-stellar mass 
loss are complicated. Here, we focus on the effects of stellar 
mass loss on a single planet and a planetesimal belt. Ob- 
servations suggest that many debris discs have their inner 
edges sculpted by planets, similar to the famous example 
of Fomalhaut ^Chiang et all 120091 ). Therefore, we consider 
a planetesimal belt with an interior planet, close enough to 
the belt such it truncates the inner edge, similar to Neptune 
and the Kuiper belt in our solar system. The planet domi- 
nates the dynamics of bodies at the inner edge of the disc: 
mater ial inside of the chaotic zone surrounding the planet's 
orbit (|Wisdomlll98bT) will be cleared, due to the overlap of 
mean motion resonances. As the star loses mass the size 
of the chaotic zone increases and extra material is scattered 
from the belt. Here we use N-body simulations to investigate 
the fate of this scattered material and whether the evolution 
of this simple system post-stellar mass loss can explain the 
white dwarf observations. In Sec. [2] we describe the set up 
simulations. In Sec. Owe outline results from our initial sim- 
ulations that mimic the main sequence evolution of the belt 
and set up the initial conditions, whilst Sec. U describes our 
simulations that include stellar mass loss. In Sec.[S]we com- 
pare our simulations to the white dwarf observations and 
finally in Sec. [S] we conclude and summarise. 



2 SETUP 

In order to investigate the dynamical effects of stellar mass 
loss on planetesimal belts we consider a simplistic planetary 
system architecture. Our simulations include a planet and 
a planetesimal belt orbiting a central star that undergoes 
mass loss. The main aim is to investigate the fate of the 
planetesimals in the belt after the star has lost mass. In 
particular, we consider the feasibility of scattering enough 
bodies towards the central star in order to produce the hot 
dusty discs observed around some white dwarfs. 

The si mulat ions are performed using Mercury 
IChambersI 1 19991 ) with the RADAU integrator. They 
are set up with a star of mass M*(t), a planet of mass M p i 
on a circular orbit and N mass-less test particles, in a belt 
initially outside of the planet's orbit. Mercury was altered 
such that the central star's mass changes as a function of 
time. The test particles are distributed in semi-major axis 
from the planet's semi-major axis a p i to a max = (j) 2 ^ 3 a p i, 
the 2:1 mea n motion resonance, the same outer edge a s the 
Kuiper belt (|Truiillo fc Brownll200ll ; lAllen et al.ll200ll ). 



Typically, high mass loss rates on the AGB last for 
~ 10 J yrs. In all our simulations we consider a IMq star 
that loses | of its mass, at a constant rate, over 10 J yrs, 
however the rate of mass loss and timescale, so long as they 
are long enough to be adiabatic, will not affect the simula- 
tions. Since many particles are removed on short timescales 
we model the belt expected at the onset of the AGB phase 
by first running the simulation for the main sequence life- 
time, tMS and removing any objects classed as scattered disc 
(see later), in addition to those ejected or scattered in. Test 
particles are given randomly selected initial semi-major axis, 
between a p i and amax, eccentricity, between and e max , in- 
clination, between and imax, mean anomaly, argument of 
pericentre and longitude of ascending node, between and 
2n. Although test particles in these simulations are evenly 
distributed in semi-major axis, different radial surface den- 
sity distributions can be considered by appropriate weight- 
ing of particles. Each particle is assigned a mass based on its 
initial semi-major axis and a disc of mass M to t, distributed 
between a p i and a max , with a surface density given by 



Yi(r)dr oc r a dr, 



(1) 



where E is the mass per unit area, r is the radial distance 
in the disc and a is a parameter of the models. Since ec- 
centricities are low, a particle's semi-major axis corresponds 
approximately to its radial position. 

During the simulations, on orbits approaching close to 
the planet, interact with it and are either scattered in to- 
wards the central star or out of the system. Some test par- 
ticles receive a large enough kick that they are put directly 
onto hyperbolic orbits, become unbound and are ejected. 
Others undergo a series of scattering interactions increas- 
ing their semi-major axis and/or eccentricity. Studies of 
comets being scattered by a single planet find that when 
they reach a distance of a ga i a ctic from the star, they are 
more str ongly influenced by the galactic tide than the cen- 
tral star (|Tremainell 19931 ). where 



(^galactic 



= W 4 AU- 



(M pl (M®)) 4 / 3 



(2) 



(M»(M )) 2 /3 a pl (AUy 

At this point they either become unbound or enter the 
Oort cloud. In our simulations we assume that a similar pro- 
cess occurs and thus any bodies that go outside of a ga i act ic 
are classified as 'ejected' and removed from the simulation, 
although occasionally bodies that enter the Oort cloud may 
return on long period orbits and re-enter the inner planetary 
system. 

The ultimate fate of bodies scattered into the inner sys- 
tem, depends on the planetary system architecture. As we 
are considering an arbitrary planetary system this is un- 
known. Observations of exo-planet systems so far suggest a 
diversity of architectures. The stability of an arbitrary plan- 
etary system post-main sequence is a com plicated dynamical 
question (e.g. lDebes fc Sigurdssor] (|2002l )). In this work our 
focus is on the dynamics of the planetesimal belt and not 
the inner planetary system and therefore we merely track 
particles that are scattered into the inner system, defined as 
a test particles scattered onto an orbit with a < dm, where 
din IS a parameter of the models. Our working assumption 
is that a fraction of the bodies that are scattered in will be 
scattered further times by inner planets and some fraction 
end up close enough to the white dwarf to be tidally dis- 
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rupted. In terms of forming the hot white dwarf discs, it is 
these particles that are of interest. 

Test particles that have interacted with the planet, but 
have not, as yet, been removed we classify as being in the 
'scattered disc', similar to the Kuiper belt's scattered disc. 
The scattered disc is defined as test particles with eccentric- 
ity higher than esD- Bodies that are in the 'scattered disc' 
will eventually be removed, scattered in or ejected, if the 
simulation were to run for sufficient time. This means that 
it is not necessary to run the simulations for the entire main 
sequence lifetime in order to find the conditions in the belt 
at the start of the AGB. Although this does mean that the 
impact of mass loss on bodies in the scattered disc at the 
end of the main sequence is neglected. 

This leaves us with four potential fates for test particles: 
scattered in (SI), ejected (EJ), in the scattered disc (SD) or 
left in the belt (B). In this work we track the number of test 
particles with each of these fates, for a range of simulations 
in which the different parameters of our model are varied. 
Parameters we can change are the number of particles, N, 
planet semi-major axis, a p i, planet mass, M p i, surface mass 
density defined by a, maximum eccentricity e max , maximum 
test particle inclination, i ma x, the radius inside of which 
bodies are considered to be scattered in, ai n , the eccentricity 
above which bodies are in the scattered disc, esD and the 
simulation time before mass loss, t^s- We also consider the 
time evolution on the post-main sequence. 



3 MAIN SEQUENCE EVOLUTION 

3.1 Baseline simulation 

For our simulation we consider a set-up similar to the solar 
system's. A Neptune mass planet (M p i — Mjv ep ) is placed 
on a circular orbit at 30AU (a p j), with 500 test particles in 
a belt extending in semi-major axis from 30 to 47.6AU (2:1 
resonance) . Test particles have initial maximal eccentricities 
e m ax = 0.1 and inclinations imax = 10° similar to the cold 
Kuiper belt. Each test particle was assigned a nominal mass 
after the simulation was completed based on a disc surface 
density profile Eq.[T]and a = 1.0, take n from sub-mm obser- 
vation s of proto-planetary discs (e.g., lAndrews fc Williams! 
(I007|)). Test particles are defined as scattered into the in- 
ner system if their semi-major axis is less than a,i n , taken 
to be a in = a p i — 7rn, where ru = a p ;(M p (/3M*) 1/ ' 3 is the 
Hill's radius, and 7rn is half the separation of Neptune and 
Uranus in our solar system, to the nearest number of Hill's 
radii. 

In this section we consider this baseline simulation and 
the effect of changing some of the parameters from this set. 
Unless explicitly stated all simulations have this set of pa- 
rameters. The mass that is ejected is defined as Mej, whilst 
Msi is the mass that is 'scattered in' and Msd is the mass 
that ends up in the scattered disc. The total mass that is 
scattered, M sca tt = Mej + Msi + Msd and is often quoted 
as a fraction of the disc mass, Mbeit- 

3.2 Setting up the initial conditions in the belt 

Before the star undergoes mass loss, the initial conditions 
in the belt must be set up. This is done by running the 
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Figure 1. The fraction of the total disc mass scattered in (SI) or 
ejected (EJ) as a function of main sequence simulation time, for 
the parameters of our initial baseline simulation. 

simulation and removing all test particles that are ejected, 
scattered in or end the simulation in the scattered disc. The 
timescale on which test particles are removed depends on 
planet mass; more massive planets scatter test particles on 
shorter timescales. We investigated the timescale on which 
test particles are removed, shown for several mass planets 
in Fig. [1] We found that the number of bodies removed 
falls off with time and therefore we chose to run our initial 
simulations for tuts = 10 7 yrs, as most test particles are 
removed in this time period. However, it is clear from Fig.[T] 
that the IM9 planet continues to scatter test particles for 
the entire main sequence (on the order 10 9 yrs). The initial 
conditions in the belt will therefore depend on the main 
sequence lifetime of a particular star. By only using 10 7 yrs 
to set up the simulation, 45% of the mass that would end 
up being removed by a lM® in 10 9 yrs is missed, whilst for 
a IMiVep this fraction is merely 10%. This should be taken 
into account when comparing the results (see later). 

3.3 The effect of varying the definition of 
'scattered in' or a,i n 

In terms of the formation of the hot white dwarf discs we are 
interested in the test particles that are 'scattered in', assum- 
ing that a fraction of these interact with planets in the inner 
system and are thus scattered onto star-grazing orbits. The 
fraction that are defined as 'scattered in', however, varies 
significantly with din • 

In order to investigate the sensitivity of our results to 
ai n in Fig. [2] we changed the definition of a;„ and calculated 
the fraction of the total mass that is defined as scattered 
in (SI). Only a few test particles spend time, at any point 
during the simulation, just inside of the planet, such that 
if ai n > 0.98a,,; then Msi ~ 0.04M(, e it- These test particles 
will generally go on to be ejected. Msi is approximately 
constant for definitions of <Un between 0.8 and 0.98. In this 
region test particles interact strongly with the planet and 
hence cover the whole range of semi-major axis space. If 
there are planets in the inner system, interior to 0.8ai„, the 
dynamics of the test particles will be dominated by interior 
planets. This behaviour is relatively independent of planet 
mass or semi-major axis and we expect it scale with the 
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Figure 2. The variation in the mass scattered in with changes to 
the parameter a;„ for the initial baseline simulation (outlined in 
Sec. I3.lt . with a 1M® (dotted), lMpf ep (solid) and e m ax = 0.1 
and i m ax = 10°. The dotted and solid vertical lines show a; n = 
a pl ~ 1i~H for IMq and lMjvep respectively. 

Hill's radius (th)- The assumption of am = a p i — 1th in the 
baseline simulation falls in this region, hence the conclusion 
of these simulations will be relatively insensitive to a m and if 
there is another planet at ~ — 13th then interactions with 
this planet will pull some fraction of planetesimals defined 
as 'scattered in' into the inner region. 

These simulations show a very interesting result in 
terms of the formation of the white dwarf discs. There is a 
lack of bodies scattered onto star-grazing orbits by a single 
planet, at 30 AU. In fact in these simulations no test particle 
has a semi-major axis less than 0.5dpi. Even including the 
eccentricity, no test particles has a pericentre of less than 
0.3a p i. In order for an asteroid to be tidally disrupted and 
form the observed discs, it must be scattered onto an orbit 
with pericentre less than the tidal radius (on the order of 
Rq). It may be that such a small percentage of test particles 
are scattered further in towards the star, but that our sim- 
ulations do not find them because we have not included a 
sufficiently large number of test particles. Given this caveat, 
these simulations show that a single planet is incapable of 
producing the observed discs and that an inner planetary 
system is necessary. 

3.4 Comparison to analytic prescription 

An analytic (or semi-analytic) model is useful to understand 
the physical mechanisms causing instability and allows us 
to get scaling laws that describe the behaviour over a wide 
range of parameter space. Analytically it is expected that 
test particles orbits close to the planet will become chaotic 
due to the overlap of mean motion resonances. This con- 
dition defines t he chaotic zone within which this occurs as 
<|Wisdomlll98(ih : 

SOchaos _ „ ( M p l \ 2/7 , . 

~i^~- g vm:) ' (3) 

where C = 1.3. 

Test particles on chaotic orbits will either be ejected or 
scattered into the inner planetary system. Thus, the fraction 
of the disc mass that will be removed, M ana iytic/Mbelt can 



Figure 3. The fates of test particles during the initial tpjs sim- 
ulation, as a function of their initial semi-major axis and eccen- 
tricity. The parameters take the values for our baseline simulation 
(outlined in Sec. 13. it . Black particles remain in the belt, whilst 
red particles are ejected, blue scattered into the inner system and 
green end the tj^s = 10 7 yrs in the scattered disc. The black dot- 
dashed line shows the size of the chaotic zone (Eq. [3), whilst the 
red dashed line shows the region q < a p \ + Sa c f lao3 . 

be calculated, for a given surface density profile, assuming 
that all test particles with initial semimajor axes less than 
dchaotic = a P i + Sa chaos are removed, 

/■2tt fa p i+da ehaos 
Manalytic = / / Y}(r)rdrd8 

JO J a pl 

= KC ™v\^i) V7 ( 4 ) 

for the index in the surface density profile, a = 1 and 
K = Mbeit/"Ka p i(2 2 l s — 1), for a disc with outer edge, a max = 

This can be compared to our N-body simulations, where 
we show that most test particles with semi-major axes less 
than a c h a otic are removed, but that some test particles with 
a. > a c haotic are also removed. Fig. [jj] shows initial semi- 
major axis and eccentricities of all test particles in the base- 
line simulation, with those test particles that are scattered 
by a lMjvep planet highlighted. The higher the initial ec- 
centricity of the test particles the higher the number of test 
particles outside of the chaotic zone that are removed. This 
also applies to inclination. No test particles were ejected in 
the 10 7 yrs of this simulation as this timescale is too long 
for a Neptune mass planet to increase a test particles semi- 
major axis to greater than a ga ; actic = 15, 000AU. 

The formulation for the chaotic zone ( Eq . [3l was devel- 
oped f or bodies on circular orbits. Althoug h Guill en fc Faberl 
(|200fJ ) showed that the same formalism applies for eccen- 
tric planets, when all test particles have the forced (or 
the planet's) eccentricity, the behaviour for test particles 
with high free eccentricities (and inclinations) is differ- 
ent. There are a few sets of simulations that show that 
the c haotic zone is larg e r for eccentric or inclined bodies 
(e.g., IVeras fc Armitagel (|2004l )). but there is no analytic 
prescription. Although the formalism for the chaotic zone 
(Eq. [3J) was developed specifically in terms of semi-major 
axis, Fig. [3] shows that most of the structure in eccentric- 
ity can be described by the scattering of test particles with 
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Figure 4. The fraction of orbits that arc chaotic are a function of 
initial semi-major axis and eccentricity for orbits around a lMj\r ep 
planet on a circular orbit at 30AU. 100 randomly selected orbits 
were calculated for each grid point, using an encounter map. The 
greyscale represents the number of orbits that become chaotic, 
black representing no chaotic orbits, white all 100 orbits were 
chaotic. The upper left hand corner of the plot has not been 
calculated as the encounter map is not valid in this regime. The 
black dot-dashed line shows the width of the chaotic zone (Eq. [3jl, 
whilst the red dashed line shows the region where test particles 
have pericentres q < a p i + 8a cnaos . 

pericentres closer to the planet than the chaotic zone i.e. 
q < cipi + 8a c haos- This may be because such particles have 
strong close encounters with the planet, or because of the in- 
creased resonance width of eccentric orbits. However, it also 
seems that the full picture is somewhat more complicated, 
especially for highly eccentric orbits. 

In order to reproduce the dependence of the chaotic 
zone width on particle eccentricity, we investigated this 
behaviour a nalytically using an en c ounter map (using t he 
formalism of iHenon fc Petit! (|l986l ); iDuncan et all (|l989h V 
The encounter map treats the particles as orbiting on un- 
perturbed Keplerian orbits, except at conjunction with the 
planet where they receive impulsive perturb ations. The ef- 
fect o f the perturbations was determined by IDuncan et alj 
(|l989l ) analytically, and the long-term evolution of an or- 
bit can be determined by repeated application of the ana- 
lytical formula. In our investigations, the map was iterated 
for 1,000 synodic periods, corresponding to ~ lMyr. Orbits 
were classified as chaotic or regular based on the Fast Fourier 
Transform of the eccentricity. Regular orbits have smooth 
Fourier Transforms with a few well-defined peaks at sev- 
eral frequencies. Chaotic orbits have power distributed over 
a range of frequencies but with large fluctuations in power 
between closely separated frequencies. To quantify the fluc- 
tuations, we defined P = log 10 |e| where e is the FFT of the 
eccentricity. The fluctuations are then quantified by: 

N= (ElL 2 (P,-P l -i) 2 /(n-l)) 1/2 (5) 

where n is the number of elements in the FFT array. Visual 
inspection of sample trajectories showed that chaotic orbits 
typically have N > 0.2, so this was adopted as the criterion 
for classifying orbits as chaotic. 

Fig. U shows the percentage of orbits that become 
chaotic for a range of initial eccentricities and semi-major 
axis, similar to Fig. [3] For each point on the grid 100 orbits, 
with random initial longitudes of pericentre and mean lon- 
gitudes, were followed and the black-white scale indicates 
the percentage of these that are classified as chaotic; black 
means that 0/100 are chaotic, whilst white means that 100% 
are chaotic. This plot reproduces the size of the chaotic zone 
for zero eccentricity particles, i.e. Eq.[3] although the factor 



C is different from the C = 1.3 given in IWisdoml (|l980T ). 
However, Fig. [4] shows that a greater fraction of orbits with 
initially higher eccentricity become chaotic for semi-major 
axis larger than the chaotic zone. This fits with the be- 
haviour observed in our N-body integrations (see Fig. [3}, 
the higher the initial eccentricities or inclinations, the more 
particles that are scattered (or end up on chaotic orbits), and 
can also be approximated by the condition q < a p i+5a c haos- 

3.5 Results 

We investigated how the mass scattered by the planet on 
the main sequence, M sca tt, varies as a function of planet 
mass. Our results are shown in Fig. [5] shows the change in 
the total mass scattered, M scat t (as a fraction of the total 
disc mass, Mbeit) and a function of planet mass, for the main 
sequence simulation (i.e. without mass loss). As anticipated, 
there is an increase in M sca tt/Mbeit with planet mass, for all 
simulations. 

A numerical simulation with e m ax = 0.0 and imax = 0.0 
is compared to an analytic prescription, assuming that all 
test particles inside of the chaotic zone are scattered (see 
Eq. |4|. The numerical simulations scatter more mass than 
predicted analytically. We consider this is due in part to the 
chaoti c zone being larger than predicted by Eq. [3J for ex- 
ample IChiang et al.l (|2009l ) find that C ~ 2.0, rather than 
C ~ 1.3, although their simulations are for a mildly ec- 
centric planet. We also find a steeper dependence in M sca tt 
with planet mass than predicted in Eq. [4] We consider that 
this is because the simulations were run for the same time 
(tMs), despite the shorter scattering timescales for higher 
mass planets. If the simulations were run for longer then 
the lower mass planets would scatter a higher fraction of 
the disc mass. 

The baseline numerical simulation, with e m ax = 0.1 and 
imax — 10.0 is compared to a prescription calculated us- 
ing the encounter map. This assumes that particles in the 
belt are randomly distributed in eccentricity and semi-major 
axis. 100 particles are placed at every grid point and using 
the calculation shown in Fig. |3J the fraction of the total 
disc mass in bodies not on chaotic orbits is calculated, for 
the given surface density profile. M aca tt calculated using the 
encounter map prescription is somewhat lower than the nu- 
merical simulations, presumably because the chaotic zone 
increases in size with inclination as well as eccentricity, but 
inclined particles were not included in the encounter map. 
We conclude that lower mass planets scatter a smaller frac- 
tion of the disc mass due to the longer scattering timescales, 
as discussed above. 

Analytically it is anticipated that these results are inde- 
pendent of dpi, which was verified by simulations, however 
changing the surface density profile changes M sca tt by up to 
15%, for 0.8 < a < 2.0. 

These simulations show that there is a definite increase 
in M SC att with initial eccentricity and inclination. Of the 
values of initial eccentricity and inclination tested here, the 
numerical simulation with e m ax = 0.2 and i m ax = 20° scat- 
ters the most test particles. This implies that an analytic 
formulation for the chaotic zone for inclined or eccentric 
particles should be larger than that given in Eq. [3] Thus, 
the expected disc structure at the end of the main sequence 
evolution will be the initial disc, minus material originally 
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Figure 5. The mass removed due to scattering by a planet as a function of planet mass. The numerical simulation with e max = 0.1 
and imax = 10.0° (dotted) should be compared to the encounter map calculation (dot dashed), whilst the numerical simulation with 
emax = and i max = 0° should be compared to the analytic prescription in which all material within the chaotic zone is scattered 
(solid). A further numerical simulation with e ma x = 0.2 and imax = 20.0° is shown. The left-hand plot is for the initial 10 7 yrs simulation 
without mass loss, whilst the right-hand plot is post-stellar mass loss and lGyr of further evolution. 



inside of the chaotic zone, that was scattered during the 
initial simulation. 



4 POST-MAIN SEQUENCE EVOLUTION 
4.1 Analytic formulation 

Once the orbital distribution at the end of the main se- 
quence has been determined using the simulations of Sec.[3l 
we then studied evolution beyond the main sequence, in- 
cluding stellar mass loss. The star loses mass on timescales 
that are long compared to the orbital timescales and thus 
this should be an adiabatic process. Indeed this is seen to 
be the case for all test particles not scattered by the planet 
and the planet itself. As the star's mass decreases by a factor 
of 3, their orbital radii increase by the same factor, whilst 
their eccentricities and inclinations remain constant. This 
would happen for all particles and the planet itself; however 
as the stellar mass decreases and the ratio of the planet's 
mass to the stellar mass increases, the zone of influence of 
the planet increases. For these simulations, where the stel- 
lar mass is decreased by a factor of 3, the size of the chaotic 
zone increases by a factor of 3 2,/7 (see Eq. [3]). Analytically 
a prediction for the amount of mass scattered can be found 
by assuming that all test particles inside of the chaotic zone 
post-mass loss, but outside of its smaller pre-mass loss value, 
are scattered, given by: 

/•27T r[a pl {0)+-i 2 ' 7 oa ahaos {0)} 

Manaiytic = T,(r)rdrd0 

JO Jla pl (0) + 5a chaos (0)} 

= A'(3 2/7 -l)5a chaos (0), (6) 

where 8a c haos(0) is the initial size of the chaotic zone (Eq.[3]), 
a p i (0) is the planet's initial semi-major axis, a is taken as 1.0 
and K = M oe i t /Tva p i(0)(2 2 ^ i — 1) is a constant determined 
from the initial belt mass. 

The right-hand panel of Fig. [5] compares the amount 
of mass scattered, following mass loss and lGyr of further 
evolution, found in the numerical simulations, to the ana- 
lytic increase in the size of the chaotic zone. The numer- 



ical simulations show approximately the same dependence 
with planet mass as the analytic prescription. The simula- 
tion with e m ax = 0.0 and i max — 0.0 is closest to the ana- 
lytic prescription, whilst as anticipated the simulations with 
higher initial eccentricities and inclinations scatter more test 
particles. 

The main cause of scattering post stellar mass loss in 
these simulations is the increase in the extent of the chaotic 
region close to the planet. This extent can be estimated an- 
alytically from Eq. [3l giving Eq. [6] however this underes- 
timates M 3C att by a factor of a few if the belt is initially 
dynamically hot. 

4.2 Scattered in or ejected? 

We investigated the fate of scattered bodies as a function 
of planet mass. This is shown in the left hand panel of 
Fig. [5] The analytic formulation does not give the ultimate 
fate of bodies, so it is necessary to use N-body simulations. 
The majority of the mass that is scattered by the planet 
ends up in the inner planetary system, according to our def- 
inition of 'scattered in'. The fraction of the mass that is 
ejected increases with planet mass. This is because higher 
mass planets give test particles a much larger kick per en- 
counter and thus are more likely to scatter bodies out of 
the system after fewer encounters. Multiple encounters are 
required however to raise the test particle's eccentricity high 
enough that it is ejected. Assuming that the Tisserand pa- 
rameter is conserved, for a test particle encountering the 
planet with a semi-major axis equal to the planet's and zero 
inclination, the condition e > 0.405 must be satisfied for 
it to be ejected. No test particles in our simulations start 
with such a high eccentricity and hence several encounters, 
each of which increase the eccentricity are required before a 
particle is ejected. 

Lower mass planets, on the other hand, tend to scatter 
test particles many many times before they get ejected. This 
leads to an increased chance of a test particle having a < a, n 
at some point before it is ejected and hence in our definition 
being scattered in. There is a general trend that higher mass 
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Figure 8. A histogram of the predicted accretion rates from our 
baseline simulation, taken from the population of discs in Fig. [7] 
Discs around white dwarfs with cooling ages between lOOMyr and 
300Myr, and less than lOGyr are shown. 



planets scatter test particles on shorter timescales and there- 
fore clear any bodies from the scattered disc before the end 
of the simulation, whereas lower mass planets end the sim- 
ulations with a higher mass in the scattered disc. However, 
the lowest mass planets scatter test particles more weakly 
and thus fewer test particles end up with e > esD = 0.24, 
our definition of the scattered disc and Msd is low for the 
lowest mass planets. 

Interestingly although M sca u is very dependent on 
planet mass, the dependence with planet mass on Msi is 
weak, at least for e max = 0.1 and i ma x = 10°; see the 
right hand panel of Fig. [6] This is because the decrease in 
Msi /Mscatt with planet mass counteracts the increase in 
M 3ca tt- Msi is highest for discs with initially high eccentric- 
ities and inclinations. This is because M aca tt is also higher 
since the size of the chaotic zone increases with eccentric- 
ity/inclination (see explanation above). Further simulations 
also show that Msi is also independent of the planet's semi- 
major axis. 



5 THE RELATIONSHIP BETWEEN THESE 
SIMULATIONS AND OBSERVATIONS OF 
METAL RICH WHITE DWARFS. 

In relation to the observations of metal polluted white 
dwarfs and/or white dwarfs with dusty discs, we are in- 
terested in bodies that are scattered into the inner plan- 
etary system. We assume that some fraction, /td, of the 
material scattered into the inner planetary system is fur- 
ther scattered by interior planets and ends up on an orbit 
that approaches close enough to the star for the body to 
be tidally disrupted. It is these bodies that potentially form 
the observed discs and accrete onto the star. This fraction 
is highly dependent on the inner planetary system archi- 
tecture and will vary between individual planetary systems. 
Observations of exo-planet systems so far suggest a diversity 
of architectures. The stability of an arbitrary planetary sys- 
tem post-main sequence is a com plicated dynamical question 
(e.g. lDebes fc Sigurdssonl fcOQj )) and it is therefore beyond 
the scope of the current work to investigate in detail scat- 
tering by the inner planetary system. 

/td can be a pproximated for our s olar s ystem from pre- 
vious simulations. ILevison fc Duncan! (|l997T I find that ~1% 
of visible comets end up on sun-grazing orbits. The parti- 
cles in their simulations that leave the Kuiper belt to become 
visible comets correspond approximatel y to our definition of 
'scattered in'. lLevison fc Duncan! (|l994T l , on the other hand, 
investigate the fate of known visible comets in our solar sys- 
tem and find that ~6% end up on sungrazing orbits. The 
discrepancy between the two is l ikely to be due to the inclu - 
sion of the terrestial planets in ILevison fc Duncan! (| 1994h . 
Therefore, we adopt a fraction /td = 0.06 of Msi that ends 
up on sun-grazing orbits and is therefore tidally disrupted. 
Of course, /td should be calculated for each individual plan- 
etary system and could vary between and 1. 

Not all of the material that is scattered in close enough 
to the star to be disrupted will end up in a disc, or be ac- 
creted onto the star. The formation of a such disc has not 
been modeled in detail at present, but here we assume that 
the disruption is relatively inefficient and only a small frac- 
tion, face ~ 0.1, of the mass that reaches ~ Rq, ends up 
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Figure 7. Left panekThe accretion rate calculated using Eq.[7]and Eq. [8] for our simulations, with M p [ = lMj (crosses and solid line), 
lAfjvep (as terisks and dotted line ) and IMj (triangles and dashed line). A range of belt masses arc calculated using the population 
models of jBonsor fc Wvattll201oh , but only the median value is plotted here, with upper and lower limits for the lMjv ep case shown 
by the arrows. Each test particle that is scattered in is plotted with a discrete value of M. A straight line is then fitted to these data 
points. Right panel: The accretion rates for a population of discs with randomly selected initial belt mass, radius, cooling age and stellar 
properties, using the smoothed fit to the stochastic accretion process, as d etermined in the le ft-hand panel, but for each individual disc 
mass. These are compared to observed heavy element accretion rates from lFarihi et al ] l|2009t ) (red asterisks). 



accreting onto the star. Hence, the mass that will be ac- 
creted onto the star is given by: 

Mace. ~ face X f T D X f S I X M be lt, (7) 

where fsi = Msi /Mbelt is the fraction of the initial belt 
mass defined as 'scattered in'. Msi can either be the total 
mass that is scattered in, and then M acc is the total mass 
that is accreted over the white dwarf lifetime, or alterna- 
tively the mass scattered in within a time interval dt, in 
which case M acc is the mass accreted in the time interval 
dt. 

Spitzer near-infrared observations of white dwarfs are 
used to determine dust masses for the observed discs. Since 
discs are opaque, this is a minimum disc mass and it is un- 
clear how it relates to the total disc mass or the mass that 
must be disrupted into order to produce such an observation. 
We, therefore, chose to compare the results of our simula- 
tions to the heavy element accretion rates calculated from 
observed abundances of metals in the white dwarf atmo- 
sphere. These are calculated from observed Ca abundances, 
an assumption that the accretion is in steady state and that 
the abundance of Ca in the accreting material is approxi- 
mately solar. 

Assuming that mass must be supplied to the disc at 
least at the rate at which it accretes onto star, the results of 
these simulations can be interpreted in terms of the obser- 
vations. The rate at which mass is scattered inwards onto 
star-grazing orbits, or the predicted accretion rate, is given 
by: 

^~~^f> (8) 

where At is the time interval over which a mass M acc is 
scattered. This assumes that the accretion is a continuous 
process and that the accretion rate is determined by the 
scattering rate rather than viscous timescales in the close-in 
disc. These accretion rates could, however, be considered a 
minimum for the rate at which material must be supplied to 
the disc in order to reproduce the observed heavy element 



accretion rates onto the star. If the pollution is produced by 
the disruption of a large individual body, as suggeste d by, 
amongst others, [jura et all (|2009t ); [Debes et all (|201lf) then 
a lower scattering rate than predicted by these simulations 
is required. 

Using this formulation and these assumptions, we cal- 
culated the accretion rate from each individual scattering 
event. The timescale for scattering, dt in Eq.|8]is calculated 
as the mean of the time between the current scattering event 
and those immediately preceding and following it. Proper- 
ties of the disc were selected randomly from the main se- 
quence population of debris discs around A stars, and the 
collisionally evolved mass at the end o f the main sequence 
dete rmined, according to th e models of IWvatt et al.l (|2007l ) 
and iBonsor fc WvattJ (|201dl ). The mass left in the disc af- 
ter our initial simulations was equated with the collisionally 
evolved mass at the end of the main sequence. Collisional 
evolu tion of the disc mass in the white dwarf phase is negli- 
gible l|Bonsor fc Wvadl201Ch . The left-hand panel of Fig. [7] 
shows these accretion rates as a function of time, for a belt 
truncated by a 1M® (crosses), lMjv ep (asterisks) and IMj 
(triangles) planet. Only the belt with the median mass at 
the end of the main sequence is displayed, whilst the arrows 
show the highest and lowest mass belt in the population, for 
the lMjvep case. 

As anticipated, early in the white dwarf phase many 
particles are scattered, whilst at later times, the number of 
particles scattered as a function of time, and thus the ac- 
cretion rate, decreases. This happens slightly more slowly 
for the lower mass planets, since scattering times decrease 
with increasing planet mass. The difference between the dif- 
ferent planet masses, however, is small, compared to the 
range of accretion rates for different initial belt properties, 
or the other assumptions that went into this plot. In order 
to convert this stochastic process into a smooth decrease 
with time, a straight line was fitted to the data for each belt 
mass. These are shown for the belt with median mass by 
the solid (1M©), dotted (lM Nep ) and dashed (IMj) lines. 
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For a lA/jVep planet the slope of this line is —1.1 ± 0.04. 
Observations of metal polluted white dwarfs also show a de- 
creas e in accretion rate with cooling age of the white dwarf 
(e.g. dFarihi et alJI^OOgh - ). Using a sample of 62 white dwarfs 
from lFarihi et aL f (|2009l . l2010ij ) we found that the decrease 
in log (t coo i) of log (M) can be fitted with a straight line of 
slope of —1.3 ± 0.23. This compares well with our simula- 
tions. 

In order to compare to the observations in more detail, 
a model population was calculated. Each star in the popu- 
lation is assumed to have a lMjvep planet, a disc with an 
initial mass and radius rand omly selected from the distribu- 
tions of IWvatt et al.l |2007l ) and a randomly selected main 
sequence lifetime from typical main sequence lifetimes for 
A stars. Ages were selected evenly distributed in log space, 
as this is consistent with the spread in ages in the observed 
sample. The smoothed formula for the decrease in accretion 
rate with time shown in the left-hand panel of Fig. [7J was 
used to calculate the accretion rates shown in the right-hand 
panel. Black crosses show our simulations, whilst heavy el- 
ement accretion rates, ca lculated from observed Ca abun- 
dances jFarihi et al.ll2009l ) are shown as red asterisks. 

Our population agrees qualitatively with these obser- 
vations, in particular, as also seen in the left-hand panel of 
Fig. there is good agreement with the decrease in accre- 
tion rate with white dwarf cooling age. There is a surpris- 
ingly good agreement with the order of magnitude and range 
of the observed accretion rates. Our models show that ac- 
cretion of material onto the star will be a continuous process 
and for the values of fro and face used in these simulations 
it is not necessary to invoke discrete events to explain the 
observations. 

However, our calculated accretion rates scale with the 
parameters /td and face, which will vary between individ- 
ual planetary systems or disruption events. Our value for 
the fraction of the disrupted material accreted onto the star 
(face = 0.1), although reasonable, may well be higher or 
lower and will vary with factors such as orbital parameters 
of the body being disrupted, its composition and strength. 
The largest uncertainty, however, is in the fraction of bod- 
ies 'scattered in' that end up on star-grazing orbits, /td- 
This will vary significantly between inner planetary system 
architectures. White dwarfs with high metal accretion rates 
that stand out from the population, such as GD362, GD40 
and HS 2253+8023, most probably have a planetary system 
that is particularly efficient at scattering bodies onto star- 
grazing orbits. Other planetary systems may be less efficient 
at scattering bodies onto star-grazing orbits, or have this ef- 
ficiency reduced post-stellar mass loss. In fact, the dynamics 
of many planetary systems will be altered post-stellar mass 
loss, potentially inner planets may be scattered such that 
they collide, are ejected or enter and clear the planetesimal 
belt. For systems where these processes are relevant, our 
simple model will no longer apply. 

To determine the importance of further dynamical pro- 
cesses and whether our values for /td and face describe 
an 'average' planetary system, we need to compare the per- 
centage of white dwarfs that are me tal polluted found by the 
observations with our simulations. IZuckerman et all (|2010T ) 
find that 19% of DB white dwarfs, with temperatures be- 
tween 13,500K and 19,500K, are metal polluted with accre- 
tion rates M > 10 8 gs _1 . For a comparable sample of DA 



white dwarfs, IZuckerman et all |2003l ) only found that 5% 
had accretion rates > 10 8 gs _1 . These differences may well 
be attributed to a differenc es in the birth environm ents of 
DA and DB white dwarfs (|Zuckerman et all |2010T |. From 
our model population in Fig. [7J we calculated a histogram 
of mass accretion rates, which are shown in Fig. [S] for two 
age samples, 100Myr< t coo i < 300Myr and t coo i < lOGyr. 
The former sample corresponds approximately to tempera- 
tures betwen 13,500K and 19,500K and 66% of this sample 
have M > 10 8 gs _1 , compared to 45% of the stars with 
100Myr< t CO oi < lOGyr. 

These figures suggest that our simulations have overes- 
timated the number of systems for which this simple model 
is applicable. The discrepancy could reflect the fraction of 
planetary systems that are destabilized post-stellar mass 
loss, since if a planet is scattered into the belt, the amount of 
material scattered may initially increase, but then decrease 
at later times as the belt is rapidly cleared. Our simula- 
tions thus suggest that either instabilities are relevant for 
many planetary systems or that most planetary systems, 
post-stellar mass loss, are significantly less efficient than our 
solar system at scattering bodies onto star-grazing orbits, 
i.e. f TD < 0.06. 

For many main sequence planetary systems our model 
is too simplistic. It no longer applies if the planet is in- 
clined or eccentric or if the dynamics of the system are 
dominated by another process, for example secular reso- 
nances or binary induced Kozai cycles. Although the ev- 
idence for planetesimal belts whose i nner edge is trun- 
cated by a planet is g ood, e.g. HR4 796 (TWyatt et al.lll999l ; 
iMoerchen et~aH l201of). HD191 089 (IChurcher et all |2010| ) , 
Fomalhaut (|Kalas et all 120051 ; IChiang et al.1 120091 1. other 
explanations have b een put forward for inne r holes, such 
as pl anet formation (jKenvon fc Bromlevl 2004 ; Smith et al.1 
120091) and interactions with gas ( Besla fc Wul 20071 ). Hence, 
the results of these simulations may not apply to all main 
sequence stars. 

Our simulations only consider a simplistic model for 
stellar mass loss. We have assumed spherical symmetry as 
the natural first assumption, however there have been sug- 
gestions tha t many systems hav e asym metric mass loss e.g 
ISokerl (|200ll ); iParriott fc Alcockl (|l998h . The stellar mass in 
the current simulations changed by a factor of 3, when in 
reality this will vary depending on the initial stellar mass, 
metallicity and so forth. This potentially changes the total 
amount of material scattered (M aC att) by up to a factor of 2, 
thus increasing the spread in the calculated accretion rates 
in Fig. [JJ 

We have also ignored other effects of stellar evolution 
that may cause a decrease in the planetesimal belt mass, 
for example stellar wind drag , YORP effect or s ublimation. 
iBonsor fc Wvattl (|2010l ) and ljura fc Xul i|2010t ) show that 
sublimation has a negligible effect on bodies of pur ely sil- 
icate or mixed composition. IBonsor fc Wvattl (|2010l ) found 
that stellar wind drag leaves between 10 -6 and 1O _1 M0 of 
material inside of the planetesimal belt at the end of the 
AGB. Th ese values will be reduced further by the resonance 
trapping (|Dong et al.l2O10l ). Nonetheless, this is significantly 
less than the total mass scattered inwards during our sim- 
ulations, between 10 -4 and 1O 2 M0, hence our simulations 
show that scattering by a planet will dominate over stellar 
wind drag. 
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To summarise, this simple model shows that if every 
star were to have a planetesimal belt truncated by a planet, 
and an inner planetary system capable of scattering bodies 
onto star-grazing orbits, this will produce the observed pol- 
lution in white dwarfs. If every system is as efficient as the 
solar system at scattering bodies onto star-grazing orbits, 
then a higher fraction of white dwarfs would be metal pol- 
luted than is found in observations. Therefore, either many 
evolved planetary systems are less efficient at scattering bod- 
ies onto star-grazing orbits, or further dynamical processes 
are important. 



6 CONCLUSIONS 

In this work we address the origin of heavy elements in metal 
polluted white dwarfs and whether accretion of asteroids 
or comets can explain these observations. We have taken a 
simple model for a Kuiper-like planetesimal belt whose in- 
ner edge is sculpted by a planet. This is typical of many 
planetary systems on the main sequence. All main sequence 
stars evolve to become giants, losing a significant propor- 
tion of their mass whilst on the asymptotic giant branch to 
end their lives as a white dwarf (or for higher mass stars 
than considered here, neutron star or black hole). We have 
used N-body simulations to investigate the effects of stellar 
evolution on this simple system, with the focus of explain- 
ing observations of metal polluted white dwarfs and white 
dwarfs with close-in dusty discs. The best models for these 
systems suggest that they are produced from asteroidal or 
cometary material that is scattered inwards due to dynam- 
ical instabilities post-stellar mass loss. 

We found that for a dynamically cold system (e m ax = 
and i m ax = 0), the amount of material scattered in the 
simulations can be calculated reasonably well using an an- 
alytic formulation, shown in Eq. [5] This assumes that pre- 
mass loss, the chaotic zone, given by Eq. [3] is cleared, whilst 
post-stellar mass loss, test particles that are inside of the in- 
creased chaotic zone are scattered. For systems with higher 
initial eccentricities and inclinations, for example e max =0.1 
and imax = 10° for the "cold" Kuiper belt, the amount of 
material scattered is higher than given by this analytic for- 
mula. The fraction of the belt mass that is scattered in- 
creases with planet mass, but is independent of planet semi- 
major axis. 

Our simulations tracked test particles that are ejected, 
scattered in, end the simulation in the scattered disc or main 
belt. Our definition of 'scattered in' included all test parti- 
cles that are scattered onto orbits with semi-major axis less 
than din — dpi — 1th- If there are interior planets, our as- 
sumption is that these bodies will interact with them and 
some fraction will be scattered onto star-grazing orbits, since 
our simulations show that a single planet is insufficient to 
scatter bodies onto star-grazing orbits. It is this fraction that 
is relevant to the white dwarf observations. We found that 
lower mass planets are more likely to scatter test particles 
into the inner planetary system, whilst higher mass planets 
give each particle a larger 'kick' with a single close encounter 
and therefore are more likely to eject them. Hence, despite 
the increase in the amount of material scattered with planet 
mass, the mass that is scattered into the inner planetary 
system is relatively independent of planet mass, given the 



caveat that the belt mass is not comparable to the planet 
mass. 

In order for our simple model to explain the white dwarf 
observations enough material must be scattered into the in- 
ner planetary system to reproduce the observations. We as- 
sume that a fraction faccfrD of the material 'scattered in' 
is accreted onto the star. Given the wealth of planetary sys- 
tem architectures found in exoplanet systems, for our calcu- 
lations we take the efficiency of the solar system at scatter- 
ing Neptune encountering bodies onto sun-colliding orbits 
(frD = 0.06), and an efficiency of the disruption process 
of face — 0.1. We assume that the initial planetesimal belt 
properties are the same as those found from observations 
of debris discs around main sequence A stars, but take into 
account the collisional evolution of disc material. Accretion 
rates are calculated using Eq. [7] and Eq. [5] for a population 
of evolved planetary systems as a function of their cooling 
age as a white dwarf (see Fig. 0). These compare well with 
the observations, reproducing the correct order of magni- 
tude, approximate range and most importantly the decrease 
in accretion rate with white dwarf age. Interestingly we find 
that that stellar mass loss can explain accretion rates even 
for old (> lGyr) white dwarfs. 

In some ways this agreement is surprising given that 
the accretion rates scale with }td and hence will vary 
significantly between individual planetary systems. In fact 
our simulations overestimate the number of highly polluted 
white dwarfs; 82% of our simulations for white dwarfs with 
tcooi < 300Myr have M > 10 8 gs~ 1 , compare d to only 19% 
of DB white dwarfs JZuckerman et ail |2010| ) or 5% of DA 
white dwarfs (|Zuckerman et al.ll2003l ). There are three fac- 
tors that could reduce the fraction of white dwarfs with high 
accretion rates calculated from our simulations. Firstly the 
efficiency of scattering bodies onto star-grazing orbits may 
be reduced by the dynamical rearrangement of the planetary 
orbits. Alternatively, altered dynamics post-stellar mass loss 
could scatter planets in such a way that they clear the plan- 
etesimal belt swiftly of material and hence accretion will not 
be observed at later times. Finally, not all main sequence 
stars will have a planetary system resembling this simple 
model, i.e. containing a planetesimal belt and interior planet 
on a circular orbit. 

This work shows that a simple model with a planetesi- 
mal belt and planet is able scatter enough material inwards 
in order to reproduce the observed metal abundances in pol- 
luted white dwarfs, even for old (>lGyr) white dwarfs. In 
fact, given the observations of debris discs and planets on 
the main sequence, this model suggests that metals should 
be observed in a higher proportion of white dwarfs than 
is found by observations. Either the solar system is particu- 
larly efficient at scattering bodies onto star-grazing orbits or 
dynamical instabilities and the rearrangement of the inner 
planetary system post-stellar mass loss is crucially impor- 
tant for many evolved planetary systems. 
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